Assignment 2

Differential Gene Expression and Preliminary ORA

1 Introduction

Bulk RNA sequencing data on 24 pre-treatment breast cancer tumours performed by Wu et al. [1] was downloaded from GEO with accession GSE176078, and associated publication.

geo_id <- "GSE176078"

This dataset was sequenced using Illumina NextSeq 500 (Homo sapiens), submitted on Apr 15 2014.

Breast cancers are clinically stratified based on:

• expression of the estrogen receptor (ER),
• expression of the progesterone receptor (PR), and
• overexpression of HER2 (or amplification of HER2 gene ERBB2)

This results in the following three clinical subtypes within this dataset:

• Luminal ie. ER+; (ER+, PR+/-)
• HER2+; (HER2+, ER+/-, PR+/-)
• Triple Negative ie. TNBC; (ER-, PR-, HER2-)

Breast cancers are also stratified on bulk transcriptomic profiling via PAM50 gene signatures [2] describing five molecular subtypes: Luminal A, Luminal B, HER2-enriched, Basal-like, and Normal-like.

The ~70-80% concordance between clinical and molecular subtypes motivated this study to improve functional understanding of breast cancer in a broader and more coherent sense.

For this study, the clinical subtype conditions classified are as follows:

• HER2+ ; (HER2+, ER-, PR+/-)
• HER2+/ER+ ; (HER2+, ER+, PR+/-)
• ER+ ; (HER2-, ER+, PR+/-)
• TNBC ; (HER2-, ER-, PR-)

The distribution of clinical subtypes across the 24 samples is shown in table 1.

Table 1: Clinical breast cancer subtype splits in our dataset

	Count	Percent
HER2+	2	8.33
HER2+/ER+	2	8.33
TNBC	8	33.33
ER+	12	50.00
Total	24	100.00

This leads to 16.6666667% of samples containing HER2 expression, and 41.6666667% of samples containing ER expression, with the overlaps of these representing 8.3333333% of all samples; 50% of those with HER2, 14.2857143% of those with ER.

The raw counts were cleaned, mapped to HUGO Gene Nomenclature Committee (HGNC) Symbols, and normalized to produce final counts which will be used for the duration of this report.

Of the original 58387 genes, we were able to map and produce 28712 unversioned, unique genes, of which filtering outliers (keeping genes present in a minimum of 12 samples) resulted in 14800 genes remaining. This results in the following plots.

The density plot shows the distribution of the cleaned, filtered, and normalized counts per million across all samples (varying colours). It displays a smooth curve as the result of our pre-processing.

Figure 1: Density of minimum 12 samples filtered and normalized CPM counts for 24 samples across four clinical subtypes of breast cancer tissue

Dispersion was calculated using edgeR to describe deviation of variance from the mean. The Biological Coefficient of Variation (BCV) is dispersion-squared, and represents the mean-variance relationship among genes.

Figure 2: Biological Coefficient of Variation (Dispersion squared) for 24 samples across four clinical subtypes of breast cancer tissue gene-wise

The normalized counts initialized above are formatted as shown by the subsection displayed in table 2.

Table 2: Format of Normalized counts dataset as described above

	CID3586	CID3921	CID3941	CID3948	CID3963
WASH7P	17.5790	1.0450	14.7274	0.0000	28.5314
CICP27	1.9276	3.6553	2.2895	0.0000	2.1353
MTND2P28	175.8840	174.5844	168.6025	354.9525	450.4684
MTATP6P1	413.8903	589.8450	462.3509	1053.1526	404.6456
FAM87B	2.6751	1.9207	2.7736	1.5650	2.4742
LINC01128	26.6312	35.3648	17.6844	20.9555	21.8894

2 Differential Gene Expression

Given our normalized expression dataset contains only grouped classification via clinical subtypes, these are the obvious choice for factors. But, due to the nature of these classifiers, the names themselves represent binary pairing of HER2 and ER expression and absence (in a clinical setting). This provides some flexibility for analyzing differential expression since we really have two protein expressions, with their binary pairing as the clinical subtype classification.

The Multi-Dimensional Scaling (MDS) plot was produced using limma on the normalized data, and it shows the distance between samples factored by the four distinct clinical subtypes.

In the given MDS plot, notice that blue and red are opposites (+/- and -/+), purple is the intermediary between blue and red (+/+), whereas green represents TNBC: double negative (-/-). Although progesterone is a marker for breast cancer, its expression is not quantifiably present within the dataset, and so is not included. However, TNBC implies PR is not overly expressed as well.

In terms of these colours in relation to each other, you can see some clustering appearing, although nothing is confidently distinct.

2.1 Calculate p-values for each gene.

I redefine the model to include the clinical subtypes of breast cancer as factors, both as four separate classifiers, and the binary pairing of over/under expression of HER2 and ER.

# aggregate model design
model_design <- model.matrix(~ types_df$subtype)

# split model design (binary pairs)
splt_model_design <- model.matrix(~ types_df$her2 + types_df$er)

Given both design choices, I used edgeR’s Quasi likelihood model to create quasi-likelihood fits.

# normalized counts grouped by clinical subtype
d <- DGEList(counts=norm_matrix, 
             group=types_df$subtype)

# estimate dispersion on subtype model design
d_ <- estimateDisp(d, model_design)
qlfit <- glmQLFit(d_, model_design)

# in parallel, estimate dispersion on binary pairing model design
d_splt <- estimateDisp(d, splt_model_design)
qlfit_splt <- glmQLFit(d_splt, splt_model_design)

Now, we can use glmQLFTest to calculate the differential expression according to this model.

# fit subtype model design
qlf.subtypes <- glmQLFTest(qlfit)

# fit split model design
qlf.her2p <- glmQLFTest(qlfit_splt, 
                        coef='types_df$erTRUE')
qlf.erp <- glmQLFTest(qlfit_splt, 
                      coef='types_df$her2TRUE')

At this point, we analyze the calculated p-values and aggregate our top hits for each model.

# top tags for aggregate model design
top_hits <- topTags(qlf.subtypes, 
                    sort.by = "PValue",
                    n = nrow(normalized_counts))

# tog tags for split model design
top_her2p <- topTags(qlf.her2p,
                     sort.by = "PValue",
                     n = nrow(normalized_counts))
top_erp <- topTags(qlf.erp,
                   sort.by = "PValue",
                   n = nrow(normalized_counts))

2.1.1 How many genes were significantly differentially expressed?

For our three model designs (one aggregate, and two binary splits), the following number of genes were significantly differentially expressed before correction.

Table 3: Number of P-Values that are significantly differentially expressed for the three different model designs

	Significantly Differentially Expressed Genes
Aggregate across subtypes	1462
HER2 expression	1240
ER expression	779

2.1.2 What thresholds did you use and why?

A threshold of 0.05 for p-values is very standard. It implies that statistically we will only achieve this outcome 5% of the time assuming the null hypothesis is true (our null hypothesis being that there is no differential expression between the clinical subtypes in these genes).

2.2 Multiple hypothesis testing

edgeR’s analysis provides a FDR value, which implies Benjamini-Hochberg, however I wanted to independently verify.

corrected_top_hits <- p.adjust(top_hits$table$PValue, 
                               method = "BH", 
                               n = nrow(normalized_counts))
length(which(corrected_top_hits < 0.05))

[1] 3

So, using this correction method, only 3 genes pass correction for significant differential expression.

2.2.1 Which method did you use? and why?

As described, the False Discovery Rates provided by our analysis match my independent understanding of the Benjamini-Hochberg method.

This method is very applicable to our use-case since it is less stringent than Bonferroni when it comes to having a large number of hypotheses. It is best practice to use Benjamini-Hochberg first to see which of your genes pass correction, and if you need a more stringent corrective method, you can consider using Bonferroni. However, Bonferroni is said to be more useful when you are testing a smaller number of hypotheses which is generally not the case for RNA sequencing data.

2.2.2 How many genes passed correction?

For our two model designs, the following number of genes passed correction for each classifier.

Table 4: Number of P-Values that are significantly differentially expressed for the three different model designs

	Significantly Differentially Expressed Genes	Significantly Differentially Expressed Genes after correction
Aggregate across subtypes	1462	3
HER2 expression	1240	0
ER expression	779	11

There is clearly a large reduction in number of genes significantly differentially expressed.

For the aggregate model of clinical subtypes, only 3 genes passed correction for false discovery. No genes passed correction for HER2 expression association, whereas 11 genes passed correction for ER expression association.

The genes and their associated values are shown in the tables below.

2.3 Amount of differentially expressed genes

2.3.1 MA Plot

An MA Plot demonstrates the log-fold change in contrast to the log-concentration for our counts. I use edgeR’s MA plot functionality.

2.3.2 Volcano Plot

The volcano plot provides another visual for the spread of our differentially expressed genes and their significance, both before and after correction.

The green points are those with significant p-values before correction, with the red points being significant genes after correction.

2.4 Visualizing top hits

2.4.1 Heatmap

Aggregating our top hits across clinical subtypes into a heatmap, we can see the expression visually as shown, along with a heatmap annotation corresponding to each binary pairing represented by the subtype classifier.

Specifically for ER overexpression or lack thereof, I added an additional heatmap.

2.4.2 Do conditions cluster together? Explain.

In figure 9, there is very visibly a cluster demonstrated between heat intensity and ER overexpression or lack thereof with the genes primarily in the top half being generally highly expressed in lack of ER, and lowly expressed in overexpression of ER. Similarly, the lower half of genes primarily follow the opposite trend; lower expression on ER- and higher expression on ER+.

In figure 10, it is more difficult to see the cluster patterns. Sample 1 has very intense RNA expressions in the top half of genes, but its matching HER2+/ER- subtype lacks this at all. This implies the association goes in a broader scale than lack of ER for these genes. Despite 11 genes being associated with overexpression of ER, it does not appear to show a visible cluster on the heatmap like the aggregate design in figure 9 does pretty well.

3 Thresholded Over-Representation Analysis

With our significantly up-regulated, and down-regulated gene-sets, we will run a thresholded gene-set enrichment analysis.

We have 807 up-regulated genes, and 655 down-regulated genes for our aggregate model design across all clinical subtypes.

I query our genesets using gProfiler’s GOSt.

# up and down regulated thresholded genes
all <- gost_res(unique(c(upreg, downreg)))
# only up regulated thresholded genes
up <- gost_res(upreg)
# only down regulated thresholded genes
down <- gost_res(downreg)

Table 7: Head of thresholded up and down regulated genes from g:GOSt

p_value	term_id	term_name
0.018483	GO:0046032	ADP catabolic process
0.018483	GO:0046031	ADP metabolic process
0.018483	GO:0006096	glycolytic process
0.018483	GO:0009057	macromolecule catabolic process
0.018483	GO:0009137	purine nucleoside diphosphate catabolic process

Table 8: Head of thresholded up-regulated genes from g:GOSt

p_value	term_id	term_name
0.0151737	GO:0022613	ribonucleoprotein complex biogenesis
0.0151737	GO:0009179	purine ribonucleoside diphosphate metabolic process
0.0151737	GO:0009181	purine ribonucleoside diphosphate catabolic process
0.0151737	GO:0009056	catabolic process
0.0151737	GO:0046032	ADP catabolic process

Table 9: Head of thresholded down-regulated genes from g:GOSt

p_value	term_id	term_name
0.0666454	GO:0031000	response to caffeine
0.0666454	GO:0071871	response to epinephrine
0.0666454	GO:0071872	cellular response to epinephrine stimulus
0.0666454	GO:0032989	cellular anatomical entity morphogenesis
0.0666454	GO:0055006	cardiac cell development

3.1 Which method did you choose and why?

I chose to use g:Profiler with R since it was well documented by Ruth Isserlin [3] and I had familiarity with it through the journal entry.

g:Profiler provides an ORA (over-representation analysis) method, and this is exactly the method I want to perform on my un-thresholded and thresholded data.

3.2 Which annotation data did you use and why?

I chose the following annotation sources to get the most comprehensive results for our functional enrichment analysis relative to biological processes and gene ontology:

• GO:BP : Gene Ontology: Biological Processes
• REAC : Reactome
• WP : WikiPathways
• GO:CC : Gene Ontology: Cellular Components
• GO:MF : Gene Ontology: Molecular Function
• KEGG : Kyoto Encyclopedia of Genes and Genomes

GO:BP, REAC, and WP were recommended by Professor Isserlin, however I ran the ORA using these and found nothing jumped out from the resulting plots. As such, I revisited to expand the scope of annotation sources in hopes for a more informative result.

3.2.1 What version of the annotation are you using?

For g:Profiler version information, I used:

• g:Profiler version: e112_eg59_p19_25aa4782
• biomaRt version: 112

For the annotation sources selected in my query, I used:

• GO:BP version annotations: BioMart classes: releases/2024-10-27
• REAC version annotations: BioMart classes: 2025-2-3
• WP version 20250110
• GO:CC version annotations: BioMart classes: releases/2024-10-27
• GO:MF version annotations: BioMart classes: releases/2024-10-27
• KEGG version KEGG FTP Release 2024-01-22

3.3 How many gene-sets were returned with what thresholds?

Table X: Gene-set sizes for different maximum thresholds

	Max 100	Max 250	Max 1000	Max 100000
Up & Down	991	1646	2160	2380
Up	558	1123	1610	1829
Down	126	324	570	643

Table X: Gene-set sizes for different maximum thresholds with minimum intersection of 10

	Max 100	Max 250	Max 1000	Max 100000
Up & Down	119	514	942	1132
Up	49	314	698	886
Down	6	49	152	193

After thresholding with maximum value of 100 and minimum intersection of 10, we can see the terms remain similar.

Table X: Head of thresholded up and down-regulated genes from ORA

	term_id	term_name	p_value	significant
1	GO:0046032	ADP catabolic process	0.018483	TRUE
2	GO:0046031	ADP metabolic process	0.018483	TRUE
3	GO:0006096	glycolytic process	0.018483	TRUE
5	GO:0009137	purine nucleoside diphosphate catabolic process	0.018483	TRUE
6	GO:0019364	pyridine nucleotide catabolic process	0.018483	TRUE

Table X: Head of thresholded up-regulated genes from ORA

	term_id	term_name	p_value	significant
2	GO:0009179	purine ribonucleoside diphosphate metabolic process	0.0151737	TRUE
3	GO:0009181	purine ribonucleoside diphosphate catabolic process	0.0151737	TRUE
5	GO:0046032	ADP catabolic process	0.0151737	TRUE
7	GO:0046031	ADP metabolic process	0.0151737	TRUE
8	GO:0006096	glycolytic process	0.0151737	TRUE

Table X: Head of thresholded down-regulated genes from ORA

	term_id	term_name	p_value	significant
4	GO:0032989	cellular anatomical entity morphogenesis	0.0666454	FALSE
7	GO:0010927	cellular component assembly involved in morphogenesis	0.0666454	FALSE
48	GO:0035051	cardiocyte differentiation	0.0676004	FALSE
135	GO:0010469	regulation of signaling receptor activity	0.0937527	FALSE
4333	GO:0005796	Golgi lumen	0.0532868	FALSE

As such, we find that up-regulated genes have some significant gene-sets after thresholding, whereas our down-regulated struggles to maintain significant values.

Our combination table shows us overlap with the up-regulated genes, further supporting the significance of our up-regulation skewed results.

3.4 Comparing up-regulated and down-regulated separate vs. together

The following figures demonstrate the over-representation analysis plots for both up & down regulated, up-regulated, and down-regulated genes respectively with only a couple of points manually highlighted.

Figure 11: Over-representation analysis of Up and Down regulated significantly differentially expressed genes

Figure 12: Over-representation analysis of Up-regulated significantly differentially expressed genes

Figure 13: Over-representation analysis of Down-regulated significantly differentially expressed genes

3.4.1 How do these results compare?

Since the up-regulated gene results tend to exceed the down-regulated gene results in terms of term significance, the combined plot appears to more closely resemble the up-regulated plot.

In our case, ORA of up and down-regulated genes in combination did not have a large effect on the results, but this is because the up-regulated genes overshadow the down-regulated. If this were a different case, we could see the difference between the up- and down-regulated plots overshadowed and melted by the combination plot. As such, it is important to look at both up- and down-regulated separately.

Although, in terms of our results of this analysis, although some terms are significant according to our arbitrary threshold, none of them are strikingly so with maximum -log(adjusted p-values) only reaching up to about 3 maximum, whereas plots seen with stronger functional enrichment presence reach up much higher towards 10-16+.

4 Interpretation

4.1 Do the over-representation results support conclusions or mechanisms discussed in the original paper?

The original paper includes substantially more processing of this data in addition to the single-cell RNA sequencing data they acquired in tandem, so it is difficult to say whether this supports those results. None of my significant differentially expressed genes appear throughout the report whatsoever. The discussion and results section for the bulk RNA-sequencing is very lacklustre compared to the rest of the analysis and results found by spatial resolution of single-cell RNA sequencing. The RNA sequencing heatmap appears possibly similar to mine, although they used scaled mean interaction score and subtyped only for ER and TNBC as opposed to all clinical subtypes. I can understand this approach since the HER2 subtype presence was relatively small compared to the rest of the samples.

Conclusively, there is not much functional enrichment apparently distinguishing the classified clinical subtypes of breast cancer; although if we had to look at the most significant term(s), our analysis would draw our attention to cadherin binding, ie. GO:0045296, of which is not present in the paper’s conclusions.

4.2 Can you find evidence, ie. publications, to support some of the results that you see? How does this evidence support your results?

Continuing with the most significant GO term from our analysis–cadherin binding–there is substantial evidence that cadherins are integral to mammary gland function, and misexpression of these leads to increased motility and consequently greater likelihood of metastases [4]. E-cadherin is one of the most widely expressed tumour suppressors in breast cancer, however cadherin-binding on its own is rather vague and difficult to draw conclusions from. We can conclude, however, that cadherins are associated with breast cancer in multiple ways, and so cadherin-binding may be a distinguishable functional enrichment between clinical subtypes.

Interestingly, one of our three significantly differentially expressed genes after correction across clinical subtypes of breast cancer is ANXA8 which is associated with breast cancer [5]. Rosetti et al. suggested ANXA8 is significantly higher in ER- subtypes as opposed to ER+ subtypes [6] which warrants further analysis from our study since across our range of clinical subtypes, we noticed ANXA8 was significantly differentially expressed. So our results support the claim that ANXA8 is expressed differentially across these subtypes of breast cancer (although specifics are not provided in our analysis yet).

So ultimately, no evidence concretely supports my results, however there is some associated evidence that suggests my results are not completely incorrect.

5 Acknowledgments

This paper makes use of packages knitr [7], BiocManager [8], GEOquery [9], kableExtra [10], edgeR [11], limma [12], ComplexHeatmap [13], circlize [14], gprofiler2 [15], GSA [16], & ggplot2 [17].

6 Bibliography

1. Wu SZ, Al-Eryani G, Roden DL, Junankar S, Harvey K, Andersson A, et al. A single-cell and spatially resolved atlas of human breast cancers. Nature genetics. 2021;53:1334–47.

2. Parker JS, Mullins M, Cheang MC, Leung S, Voduc D, Vickery T, et al. Supervised risk predictor of breast cancer based on intrinsic subtypes. Journal of clinical oncology. 2009;27:1160–7.

3. Isserlin R. Run g:profiler from r.

4. Derycke LD, Bracke ME. N-cadherin in the spotlight of cell-cell adhesion, differentiation, embryogenesis, invasion and signalling. International Journal of Developmental Biology. 2004;48:463–76.

5. LifeMap Sciences Inc. MalaCards human disease database.

6. Rossetti S, Bshara W, Reiners JA, Corlazzoli F, Miller A, Sacchi N. Harnessing 3D models of mammary epithelial morphogenesis: An off the beaten path approach to identify candidate biomarkers of early stage breast cancer. Cancer letters. 2016;380:375–83.

7. Xie Y. Dynamic documents with R and knitr. 2nd edition. Boca Raton, Florida: Chapman; Hall/CRC; 2015.

8. Morgan M, Ramos M. BiocManager: Access the bioconductor project package repository. 2024.

9. Davis S. GEOquery: Get data from NCBI gene expression omnibus (GEO). 2024.

10. Zhu H. kableExtra: Construct complex table with ’kable’ and pipe syntax. 2024.

11. Chen Y, Lun AT, McCarthy DJ, Chen L, Baldoni P, Ritchie ME, et al. edgeR: Empirical analysis of digital gene expression data in r. 2025.

12. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research. 2015;43:e47.

13. Gu Z. Complex heatmap visualization. iMeta. 2022. https://doi.org/10.1002/imt2.43.

14. Gu Z, Gu L, Eils R, Schlesner M, Brors B. Circlize implements and enhances circular visualization in r. Bioinformatics. 2014;30:2811–2.

15. Kolberg L, Raudvere U, Kuzmin I, Vilo J, Peterson H. gprofiler2– an r package for gene list functional enrichment analysis and namespace conversion toolset g:profiler. F1000Research. 2020;9 (ELIXIR).

16. Efron B, Tibshirani R. GSA: Gene set analysis. 2024.

17. Wickham H, Chang W, Henry L, Pedersen TL, Takahashi K, Wilke C, et al. ggplot2: Create elegant data visualisations using the grammar of graphics. 2024.